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ABSTRACT 

We examine the thermodynamical properties of a family of partially relaxed, 
anisotropic stellar systems, derived earlier from the Boltzmann entropy under 
the assumption that a third quantity Q is conserved in addition to the total 
energy and the total number of stars. We now show that the family of models 
conforms to the paradigm of the gravothermal catastrophe, which is expected 
to occur (in the presence of adequate energy transport mechanisms) when the 
one-parameter equilibrium sequence attains sufficiently high values of the concen- 
tration parameter; these are the values for which the models are well fitted by the 
i?^/^ law. In the intermediate concentration regime the models belonging to the 
sequence exhibit significant deviations from the i?^/^ law. Curiously, in the low- 
concentration regime, the global thermodynamical temperature associated with 
the models becomes negative when the models become too anisotropic so that 
they are unstable against the radial orbit instability; this latter behavior, while 
offering a new clue to the physical interpretation of the radial orbit instability, 
is at variance with respect to the low-concentration limit of the classical case of 
the isotropic, isothermal sphere investigated by Bonnor (1956) and Lynden-Bell 
and Wood (1968). 



Subject headings: stellar dynamics — galaxies: evolution — galaxies: formation 
— galaxies: kinematics and dynamics — galaxies: structure 
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1. Introduction 

The possibility of providing a thermodynamical description of self-gravitating stellar 
systems has motivated a number of investigations in galactic dynamics, starting with the 
pioneering work of Antonov and Lynden-Beli in the 60s. After the realization that vio- 
lent relaxation is likely to lead to partially relaxed configurations in dynamical equilibrium 
(Lynden-Bell 1967), a re-examination of the problem of the isothermal sphere, studied earlier 
by Bonnor (1956) for a self-gravitating gas, led to the interesting possibihty that stellar sys- 
tems may undergo the process of gravothermal catastrophe (Lynden-Bell and Wood 1968; 
see also Antonov 1962). However, a rigorous derivation of the onset of the gravothermal 
catastrophe from a study of the Boltzmann entropy appeared to be available only for the 
case of ideal systems confined by a spherical refiecting wall. A number of convincing qual- 
itative arguments made it clear that also unbound stellar systems with finite mass, such 
as those described by the King sequence (King 1966; these spherical models have a finite 
radius, but do not require an external wall), should probably fall into the same physical 
framework and indeed the paradigm received a lot of attention, especially in the context of 
the dynamics of globular clusters (see Spitzer 1987), which are known to possess, at least 
to some extent, the desired internal coUisionality (see also Lynden-Bell and Egglcton 1980). 
[An indirect indication that the general physical picture of the gravothermal catastrophe is 
likely to be robust comes also from the proof that the behaviour of the classical gas case is 
basically independent of the assumption of spherical geometry (Lombardi and Bcrtin 2001).] 
Several investigations have aimed at producing a rigorous derivation of the gravothermal 
catastrophe for unbound stellar systems, focusing on the underlying argument that refers to 
the Poincare stability of linear series of equihbria (Katz 1978, 1979; Padmanabhan 1989), 
but the proof has always been centered on an unjustified ansatz in order to connect the 
underlying entropy S with the global temperature T = l/{dS/dEtot) (see Appendix V in the 
article by Lynden-Bell and Wood 1968; Katz 1980; Magliocchetti et al. 1998). Other investi- 
gations have explored the possibility of setting the discussion in the context of non-standard 
entropies (e.g., the Tsallis entropies; see the study of the polytropic spheres by Chavanis 
2002 and Taruya and Sakagami 2002). Note that the concept of entropy for collisionless 
systems is quite subtle (e.g., see StiaveUi and Bertin 1987 and references therein). One might 
even argue whether it is actually compatible with the long-range nature of gravity, given the 
fact that self-gravitating systems lack additivity, a key ingredient in thermodynamics. 

In the meantime, inspired by N-body simulations of collisionless collapse (van Albada 
1982), which confirmed the general picture of incomplete violent relaxation and showed that 
it can lead to systems with realistic density profiles without ad hoc tuning of the initial 
conditions, some families of models were constructed able to reproduce, for quasi-spherical 
configurations, the characteristic feature of the anisotropy profile with an inner isotropic core 
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and an outer radially biased envelope (Bertin and Stiavelli 1984; see Bertin and Stiavelli 
1993 and references therein): these families turned out to exhibit the characteristic R^^^ 
projected density profile and indeed were shown to match nicely the observed photometric 
and kinematical characteristics of bright ellipticals. In an attempt at providing a justification 
of these models (in particular, of the so-called foo models, constructed initially only from 
dynamical arguments) from statistical mechanics, two routes were pursued (Stiavelli and 
Bertin 1987). 

The first combines an explicit statement of partial relaxation, i.e. of a relaxation process 
that is expected to be inefficient in the outer regions, and the existence of a suitable weight, 
related to the orbital period, for the cells that make the relevant partition of phase space; 
it follows qualitative arguments proposed by Lynden-Bell (1967) and is physically appealing 
(see also Tremaine 1986). It was indeed shown to lead naturally to the foo models. However, 
this route is not fully satisfactory from the mathematical point of view, especially since it 
involves an explicit approximation for the orbital period that is applicable only to the low 
binding energy limit of quasi-Keplerian orbits. The second route is straightforward from 
the mathematical point of view, being based on the classical Boltzmann entropy and on the 
assumed explicit conservation of a third quantity Q, in addition to the total mass M and 
to the total energy Ef^f It was shown to lead to an analytically different family of models 
(that we may call the f^"^ models; see definition in Sect. 3 below), with quafitative properties 
very similar to those of the foo models. Those models were not studied much further and 
did not receive great attention, not only because the relevant distribution function is not as 
simple as that of the foo models, but especially because the conservation of Q could only be 
justified approximately by inspection of a number of N-body simulations, without a clear- 
cut physical justification (see Stiavelli and Bertin 1987; in contrast, the conservation of the 
additional A ■ B invariant sometimes invoked in plasma physics is rather transparent; see 
Chandrasekhar and Woltjer 1958). 

In this paper we take advantage of the simple statistical mechanics foundation of the f^"^ 
family of unbound partially relaxed stellar systems to explore the possibility of a thermody- 
namical description of stellar-dynamical models that are endowed with realistic properties. 

2. Comment on the evolution of elliptical galaxies 

Before proceeding to illustrate the results of this paper, we make a short digression in 
order to bring out the connections between the present analysis and the evolution of el- 
liptical galaxies. We start by recalling that, formally, the sequence of King (1966) models 
is one special family of solutions of the coUisionless Boltzmann equation. Yet, it is recog- 
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nized to provide a reasonable description of the current properties of globular clusters (see 
Djorgovski and Meylan 1994), within a framework where these stellar systems continually 
evolve as a result of a variety of mechanisms (among which star evaporation and disk shock- 
ing; see Vesperini 1997 and references therein) and where the paradigm of the gravothermal 
catastrophe can be applied (see Spitzer 1987). Of course, it is well known that the level of 
internal coUisionality in globular clusters is relatively high, so that the above approach is 
quite natural. 

In contrast, one might at first think of dismissing the possibility that the paradigm 
of the gravothermal catastrophe should be of interest for the study of elliptical galaxies, 
because these large stellar systems lack the desired level of coUisionality, judging from the 
estimate of the relevant star-star relaxation times. Here, following the spirit of earlier inves- 
tigations (starting with Lynden-Bell and Wood 1968), we note that real elliptical galaxies 
are actually complex systems the evolution of which goes well beyond the idealized frame- 
work of the collisionless Boltzmann equation. In other words, splitting their description into 
past (formation) processes and present (mostly collisionless equilibrium) conditions should 
be considered only as an idealization introduced in order to assess the properties that define 
their current basic state. 

In practice, elliptical galaxies arc expected to be in a state of continuous evolution, for 
which we can list several specific dynamical causes: (1) Left-over granularity of the stellar 
system itself from initial collapse. Clumps of stars are likely to continue to form and dissolve 
in phase space even after the system has reached an approximate steady state. This acts as 
internal coUisionality thus making some relaxation proceed even at current epochs. Indeed, 
numerical simulations of violent (partial) relaxation show that some evolution continues well 
after the initial collapse has taken place. (2) Drag of a system of globular clusters or other 
heavier objects towards the galaxy center. A globular cluster system or the frequent capture 
of small satellites (mini-mergers) may provide an internal heating mechanism associated with 
the process of dynamical friction by the stars on the heavier objects (Bertin at al. 2002). (3) 
Long-term action of tidal interactions of the galaxy with external objects. (4) Presence of 
gas in various phases (cold, warm, and hot). Significant cooling fiows have been observed in 
bright ellipticals. Traditionally, studies of processes of this kind focus on the dynamics of the 
cooling gas and keep the background stellar system as 'frozen'. In reality, energy and mass 
exchanges take place between the stellar system and the interstellar medium. (5) Interaction 
between the galactic nucleus and the galaxy. A number of interesting correlations have been 
found between the properties of galaxy nuclei and global properties of the hosting galaxies 
(e.g., see Pellegrini 1999). These correlations suggest that significant energy exchanges are 
taking place between the galaxy and its nucleus. Eventually, if a sufficiently concentrated 
nucleus is generated, then star-star relaxation in the central regions may also become a 
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significant cause of dynamical evolution. 

All of the above are specific mechanisms that are expected to make eUiptical galaxies 
evolve in spite of their very long typical star-star relaxation time. Most of these processes are 
hard to model and to calculate in detail. As for the evolution of other complex many-body 
systems, it is hoped that thermodynamical arguments may help us identify general trends 
characterizing such evolution. This is the basic physical scenario in which the calculations 
presented in this paper are expected to be of interest for real elliptical galaxies. 



3. Partially relcLxed, unbound, finite mass systems from the Boltzmann 

entropy 

Let us consider the standard Boltzmann entropy S — — J / In fd^x(Pv and look for func- 
tions that extremize its value under the constraint that the total energy Efo^ — (1/3) J Efd^xd^v, 
the total mass M — J fd^xd^v, and the additional quantity 

Q^JriEr-'Wv (1, 

are taken to be constant. Here the functions E and represent specific energy and specific 
angular momentum square of a single star subject to a spherically symmetric mean potential 
$(r). As shown elsewhere (Stiavelli and Bertin 1987), this extremization process leads to 
the following family of distribution functions 



= Aexp 



-aE-d 



J2 



|E|3/2 



(2) 



where a, A, and d are positive real constants. One may think of these constants as pro- 
viding two dimensional scales (for example, M and Q) and one dimensionless parameter; 
the dimensionless parameter can be taken to be 7 = ad"^^^ /{4:TtGA). In principle, v is any 
positive real number; in practice, we will focus on vahies of /y ~ 1. The f^''^ non-truncated 
models are constructed by taking this form of the distribution function for E < 0, a. van- 
ishing distribution function for E > 0, and by integrating the relevant Poisson equation 
under the condition that the potential $ be regular at the origin and behaves like —GM/r 
at large radii. This integration leads to an eigenvalue problem (see Appendix) for which a 
value of 7 is determined by the choice of the central dimensionless potential, 7 = 7(^'), with 
^ = -a$(r = 0). 
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The main point of the following analysis is the determination of the Boltzmann entropy 

S'(M, Q,^) and of the total energy Etot{M^Q^'^) along the sequence of models, i.e. as a 
function of the concentration parameter \E' defined above. These functions, at constant M 
and Q, are illustrated in Fig. 1. They have been obtained by noting that, from the definitions 
of 5 and f^^\ 

S = -M In A + 2,aEtot + dQ. (3) 

From the definitions Q = Aa-^/'^d-'^-^/''Q{^) and M = Aa'^/^d-^/" M{'<iJ) and the definition 
of 7, we can express the variables {A,a,d) in terms of the variables (M, Q, ^1/) and thus find 
that the entropy per unit mass can be written as S/M — Sq{M,Q) + a{^), where Sq is 
constant when the values of M and Q are fixed, with 

(7 = -In M~(35^7-5 + — + . (4) 
V /MM 

Here E = Ei^^) is the dimcnsionless total energy defined from Etot = Aa^^^^^d^^^'^E. From 
the identity aEtot/M = E/M and the expression of a = a(M, Q, ^) obtained previously, we 
find Etot/M = H{M, g)e(*), with: 

4 ^ 9i'-\-4 ^ 4 ^ 

e = ^5M 5irQ5^E . (5) 

The factor H{M, Q) is a constant when M and Q are taken to be constant. The quantities 
7(^), M{^), Q{^), and E{^) that enter the expression of a and e depend only on ^ and 
are evaluated numerically on the equilibrium sequence. 

This completes the derivation that allows us to draw the analogy with the classical paper 
of Lynden-Bell and Wood (1968). This step, straightforward for the f^^'> models, is by itself 
interesting and new. In fact, other attempts at applying the paradigm of the gravothermal 
catastrophe to stellar dynamical equilibrium sequences were either based on an unjustified 
ansatz for the identification of the relevant temperature (e.g., see Appendix V in the article 
by Lynden-Bell and Wood 1968; Katz 1980; Magliocchetti et al. 1998) or on the use of 
non-standard entropies (for less realistic models; Chavanis 2002). 
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4. The high concentration regime: gravothermal catastrophe 

When the f^'^^ models were constructed (Stiavelh and Bertin 1987), it was immediately 
realized that they have general properties similar to those of the foo models (Bertin and 
Stiavelli 1984); in particular, for values of z/ ~ 1, sufficiently concentrated models along the 
sequence tend to settle into a "stable" overall structure, except for the development of a more 
and more compact nucleus, as the value of increases, and are characterized by a projected 
density profile very well fitted by the i?^/^ law characteristic of the surface brightness profile 
of bright elliptical galaxies. This property is illustrated in Fig. 2. 

Now, by inspection of Fig. 1 and by analogy with the study of the isothermal sphere 
(Lynden-Bell and Wood 1968), we can identify the location at ^ 9 as the location for 
the onset of the gravothermal catastrophe. This sequence of models thus has the surprising 
result that the value of that defines the onset of the gravothermal catastrophe is precisely 
that around which the models appear to become realistic representations of bright elliptical 
galaxies. We leave to other papers (see Bertin and Stiavelli 1993) the detailed discussion of 
the issues that have to be addressed when comparison is made with the observations. 

We note that in this regime of high concentration the general properties of the gravother- 
mal catastrophe are reasonably well recovered by the use of the ansatz that the temperature 
parameter conjugate to the total energy is a, a quantity directly related to the velocity 
dispersion in the central regions. Basically, this was the ansatz made in the discussion of 
the possible occurrence of the gravothermal catastrophe for the King models or for other 
sequences of models (e.g., see Lynden-Bell and Wood 1968, Katz 1980, Maghocchetti et al. 
1998). Here we have proved that the apphcation of a rigorous derivation, which is available 
in our case, gives rise to relatively modest quantitative changes in the (£'tot, ^/T) diagram 
for values of \l/ close to and beyond the onset of the catastrophe (see Fig. 3). However, in 
Sect. 6 we will draw the attention to an interesting, qualitatively new phenomenon missed 
in the previous derivations based on the use of the a- ansatz. 

In passing, we note that in this regime of relatively high concentrations, the /^''^ models 
possess one intrinsic property that makes them more appealing than the widely studied /oo 
models. This is related to the way the models compare to the phase space properties of the 
products of collisionless collapse, as observed in N-body simulations (van Albada 1982). In 
fact, one noted unsatisfactory property of the concentrated /oo models was their excessive 
degree of isotropy with respect to the models produced in the simulations.^ Here we can 



^Merritt et al. (1989) stressed this point and thus argued that a better representation of N-body simula- 
tions would be obtained by considering the /oo family of models extended to the case of negative values of 
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easily check that the anisotropy level of the concentrated f^"^ models, while still within the 
desired (radial orbit) stability boundary and still consistent with the modest amount of radial 
anisotropy revealed by the observations, seems much closer to that resulting from N-body 
simulations of collisionless collapse; in particular, the anisotropy radius r^, defined from the 
relation a{rQ.) = 1, with a = 2 — {{vj) + {v^), is close to the half-mass radius tm (while 
for the /oo models it is about three times as large). This is illustrated in Fig. 4. In any 
case, we would like to emphasize that, while we have been clearly taking inspiration from 
simulations of collisionless collapse, our main interest is in comparing the structure of our 
models with that of observed objects rather than in providing a detailed fit to the results of 
N-body simulations. 

5. The intermediate concentration regime: the R^^* law and deviations from it 

The intermediate concentration regime (the precise point that marks the low concen- 
tration regime will be identified in the next Section) is a regime where the models appear to 
be stable, with respect not only to the gravothermal catastrophe (following the arguments 
provided earlier; but we should recall that the catastrophe is expected to require a suffi- 
ciently high level of effective coUisionality in order to take off) but also to other instabilities 
(see the discussion given by Bcrtin et al. 1994 and references therein). The relatively wide 
variations, between ^ = 3.5 and ^ = 9, in all the representative quantities that characterize 
the equilibrium models suggest that this part of the sequence could be used to model the 
weak homology of bright elliptical galaxies (see Bertin et al. 2002), much like the sequence 
of King models is able to capture observed systematic variations in the structure of globular 
clusters (see Djorgovski and Meylan 1994). 



the coefficient a multiplying the energy in the exponent. Unfortunately, their proposed solution, in terms of 
models characterized by such a peculiar phase-space structure, is imablc to reproduce both the core velocity 
distribution observed in numerical experiments and the modest amount of radial anisotropy revealed by 
observed line profiles. In addition, their proposed solution is not viable because for negative values of a the 
radial anisotropy level is so high that the models are violently unstable, on an extremely short timescale, 
with the result that their structural properties would be drastically changed by rapid evolution (Stiavelli 
and Sparke 1991; Bertin et al. 1994). 
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6. The low concentration regime: negative global temperature and radial 

orbit instability 

The low concentration regime is marked by an unexpected and significant difference 
with respect to the low concentration limit of the classical isothermal sphere (Bonnor 1956; 
Lynden-Bell and Wood 1968). In fact, while the classical case reduces to the ideal non- 
gravitating gas, to which Boyle's law apphes, for the /*^^) models the system remains self- 
gravitating, although it develops a wide core in the density distribution. A clear-cut proof 
of this difference is given by inspecting the behavior of the global temperature T, identified 
from the thermodynamical definition T = 1/ {dS/ dEtot)- While the temperature defined by 
the a-ansatz remains obviously positive definite, by definition, if we look at Fig. 1 we see 
that the global temperature T changes sign at ^ a; 3.5. This marks a drastic qualitative 
deviation from the classical studies. 

Here we note a curious coincidence of this transition value of \E' with the value around 
which the sequence is bound to change its stability properties with respect to the radial orbit 
instability (Polyachenko and Shukhman 1981). Indeed, the location where the sequence is 
expected to become unstable in this regard is precisely that defined by \I' ~ 3.5, as can be 
judged from inspection of Fig. 4; around those values of ^ the level of radial anisotropy, as 
measured by 2Kj./Kt reaches the threshold value of 1.8 — 2, known to be sufficient for the 
excitation of the instability (the precise value of 2^^/ Kt corresponding to marginal stability 
is model dependent; for some sequences the reported value is below the range suggested by 
Polyachenko and Shukhman 1981). 

These clues appear to be interesting and important, but more work is required before a 
final claim can be made that there is indeed a direct relation between the dynamical radial 
orbit instability and the fact that the system possesses a negative global temperature, as we 
found based on the simple work presented here. 

7. Conclusions 

A relatively straightforward and simple thermodynamical description of an equilibrium 
sequence constructed earlier and known to possess realistic characteristics with respect to 
bright elliptical galaxies shows that, on the one side, the paradigm of the gravothermal 
catastrophe may be adequate to explain the occurrence of realistic properties in models of 
collisionless stellar systems and, on the other side, a long-known dynamical instability might 
turn out to be interpreted in terms of a thermodynamical argument. 

Probably the main open question regarding the models discussed in this paper, partially 
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addressed in previous studies (see Stiavelli and Bertin 1987), is to what extent the quantity 
Q is actually reasonably well conserved during violent (partial) relaxation. It is likely that 
a thorough investigation of this issue may give indications that the quantity is best con- 
served only in certain ranges of v. Studies of this type, combined with other dynamical and 
thermodynamical considerations, may turn out to lead to the identification of a family of 
equilibrium models with optimal behavior with respect to statistical mechanics, with respect 
to what we know about coUisionless collapse, and with respect to the problem of providing 
a realistic representation of bright elliptical galaxies. 



We would hke to thank L. Ciotti, M. Stiavelli, and T. van Albada for interesting dis- 
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A. Numerical integrations 

The f'^"'' models are constructed by solving the Poisson equation 

44:r'^$(r) = -P(r,<i>), (Al) 
dr dr 7 

in which 7 is considered as an eigenvalue to be determined by imposing the two natural 
boundary conditions $(0) = and l>(f) ~ — M(f)/(47r7f) as f — cxd. Here the symbol" 
indicates that the quantity is suitably expressed in dimensionless form. 



We have computed the two-dimensional integral for the density with an adaptive seven- 
point scheme (Berntsen et al. 1991) in order to properly handle the presence of a peaked 
integrand for certain values of the pair (f, <l>). The Poisson equation has then been solved 
with a fourth order Runge-Kutta code by starting from f = with a seed value for 7 and 
iterating the procedure until the boundary condition at large radii is matched within a certain 
accuracy. Finally, we have proceeded to calculate the global quantities M{^), Q{^), and 
E{^!) from their definitions. 



In order to check the accuracy of the numerical integration we have performed the following 
tests: (1) The virial theorem is satisfied with accuracy of the order 10"^ or better; (2) The 
integrated mass (from its definition) and the mass derived from the asymptotic behaviour of 
the potential at large radii are the same with accuracy of the order 10^"^; (3) The expression 
for <l>(f) at large radii to two significant orders in the relevant asymptotic expansion has 
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been checked to be correct with an accuracy from 10~^ to 10~^; (4) The asymptotic analysis 
allows us to estimate the contributions to M, Q, and E external to a sphere of large radius 
i?; this has been checked to help improve the numerical determination of the relevant global 
quantities. 

We estimate that the final relative error in the quantities along the equilibrium sequence 
is of the order of some parts times 10"*^ for M and Q and some parts times 10~^ for E. The 
total energy is less sensitive to the finite radius truncation error, due to its 1/ convergence. 
The propagation of these errors leads to the error bars plotted in Fig. 1, where the entropy 
a is the more difficult to determine with good accuracy. 
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Fig. 1. — Specific entropy and total energy along the equifibrium sequence of /^^^^ models 
with v — 1 (as a function of the concentration parameter ^, at constant M and Q, and 
thus expressed by means of the functions (7(^) and defined in the text). Note that for 
^ ^ 3.5 the models are characterized by a negative temperature, because the derivatives of 
S and Etot have opposite signs. 
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Fig. 2. — Residuals fi^'^^ — jj}/'^ obtained by fitting the R^^'^ law to the projected density 
profile of /^'^^ models lor v — 1 and some values of . 
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Fig. 3. — Instability spiral of f^'^^ models with v = 1. The solid line refers to the results 
obtained with the a-ansatz (with a = Crosses represent the global 

temperature from the definition dS/dEtot] other symbols indicate estimated points for which 
the adopted numerical differentiation is less reliable. The values of ^ and e for points A and 
B with a vertical tangent remain unchanged. 
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Fig. 4. — Anisotropy along the equilibrium sequence: anisotropy radius in units of the half 
mass radius 2ra/rM and the anisotropy parameter 2Kr/KT (ratio of total kinetic energy in 
the radial direction to that in the tangential directions) of f^'^^ models with v — 1. 



